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We present a new learning algorithm for Boltzmann machines that 
contain many layers of hidden variables. Data-dependent statistics are 
estimated using a variational approximation that tends to focus on a sin- 
gle mode, and data-independent statistics are estimated using persistent 
Markov chains. The use of two quite different techniques for estimating 
the two types of statistic that enter into the gradient of the log likelihood 
makes it practical to learn Boltzmann machines with multiple hidden lay- 
ers and millions of parameters. The learning can be made more efficient 
by using a layer-by-layer pretraining phase that initializes the weights 
sensibly. The pretraining also allows the variational inference to be ini- 
tialized sensibly with a single bottom-up pass. We present results on the 
MNIST and NORB data sets showing that deep Boltzmann machines 
learn very good generative models of handwritten digits and 3D objects. 
We also show that the features discovered by deep Boltzmann machines 
are a very effective way to initialize the hidden layers of feedforward 
neural nets, which are then discriminatively fine-tuned. 


1 A Brief History of Boltzmann Machine Learning 

The original learning procedure for Boltzmann machines (see section 2) 
makes use of the fact that the gradient of the log likelihood with respect 
to a connection weight has a very simple form: it is the difference of two 
pair-wise statistics (Hinton & Sejnowski, 1983). The first statistic is data 
dependent and is the expectation that both binary stochastic units in a pair 
are on when a randomly selected training case is clamped on the "visible" 
units and the states of the "hidden" units are sampled from their posterior 
distribution. The second statistic is data independent and is the expectation 
that both units are on when the visible units are not constrained by data 
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and the states of the visible and hidden units are sampled from the joint 
distribution defined by the parameters of the model. 

Hinton and Sejnowski (1983) estimated the data-dependent statistics by 
clamping a training vector on the visible units, initializing the hidden units 
to random binary states, and using sequential Gibbs sampling of the hidden 
units (Geman & Geman, 1984) to approach the posterior distribution. They 
estimated the data-independent statistics in the same way, but with the ran- 
domly initialized visible units included in the sequential Gibbs sampling. 
Inspired by Kirkpatrick, Gelatt, and Vecchi (1983), they used simulated an- 
nealing from a high initial temperature to a final temperature of one to 
speed up convergence to the stationary distribution. They demonstrated 
that this was a feasible way of learning the weights in small networks, but 
even with the help of simulated annealing, this learning procedure was 
much too slow to be practical for learning large, multilayer Boltzmann ma- 
chines. Even for small networks, the learning rate must be very small to 
avoid an unexpected effect: the high variance in the difference of the two 
estimated statistics has a tendency to drive the parameters to regions where 
each hidden unit is almost always on or almost always off. These regions 
act as attractors because the variance in the gradient estimate is lower in 
these regions, so the weights change much more slowly. 

Neal (1992) improved the learning procedure by using persistent Markov 
chains. To estimate the data-dependent statistics, the Markov chain for each 
training case is initialized at its previous state for that training case and 
then run for just a few steps. Similarly, for the data-independent statistics, 
a number of Markov chains are run for a few steps from their previous 
states. If the weights have changed only slightly, the chains will already 
be close to their stationary distributions, and a few iterations will suffice 
to keep them close. In the limit of very small learning rates, therefore, the 
data-dependent and data-independent statistics will be almost unbiased. 
Neal did not explicitly use simulated annealing, but the persistent Markov 
chains implement it implicitly, provided that the weights have small initial 
values. Early in the learning, the chains mix rapidly because the weights are 
small. As the weights grow, the chains should remain near their stationary 
distributions in much the same way as simulated annealing should track 
the stationary distribution as the inverse temperature increases. 

Neal (1992) showed that persistent Markov chains work quite well for 
training a Boltzmann machine on a fairly small data set. For large data sets, 
however, it is much more efficient to update the weights after a small mini- 
batch of training examples, so by the time a training example is revisited, the 
weights may have changed by a lot and the stored state of the Markov chain 
for that training case may be far from equilibrium. Also, once the weights 
become large, the Markov chains used for estimating the data-independent 
statistics may have a very slow mixing rate since they typically need to 
sample from a highly multimodal distribution in which widely separated 
modes have very similar probabilities but the vast majority of the joint states 
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are extremely improbable. This suggests that the learning rates might need 
to be impractically small for the persistent chains to remain close to their 
stationary distributions with only a few state updates per weight update. 
Fortunately, the asymptotic analysis is almost completely irrelevant: there 
is a subtle reason, explained later in this section, why the learning works 
well with a learning rate that is much larger than the obvious asymptotic 
analysis would allow. 

In an attempt to reduce the time required by the sampling process, Peter- 
son and Anderson (1987) replaced Gibbs sampling with a simple mean-field 
method that approximates a stationary distribution by replacing stochastic 
binary values with deterministic real-valued probabilities. More sophisti- 
cated deterministic approximation methods were investigated by Galland 
(1991) and Kappen and Rodriguez (1998), but none of these approximations 
worked very well for learning for reasons that were not well understood at 
the time. 

It is now well known that in directed graphical models, learning typically 
works quite well when the statistics from the true posterior distribution 
that are required for exact maximum likelihood learning are replaced by 
statistics from a simpler approximating distribution, such as a simple mean- 
field distribution (Zemel, 1993; Hinton & Zemel, 1994; Neal & Hinton, 1998; 
Jordan, Ghahramani, Jaakkola, & Saul, 1999). The reason learning still works 
is that it follows the gradient of a variational bound (see section 2.2). This 
bound consists of the log probability that the model assigns to the training 
data penalized by the sum, over all training cases, of the Kullback-Leibler 
divergence between the approximating posterior and the true posterior 
over the hidden variables. Following the gradient of the bound tends to 
minimize this penalty term, thus making the true posterior of the model 
similar to the approximating distribution. 

An undirected graphical model, such as a Boltzmann machine, has an 
additional, data-independent term in the maximum likelihood gradient. 
This term is the derivative of the log partition function, and unlike the 
data-dependent term, it has a negative sign. This means that if a varia- 
tional approximation is used to estimate the data-independent statistics, 
the resulting gradient will tend to change the parameters to make the ap- 
proximation worse. This probably explains the lack of success in using 
variational approximations for learning Boltzmann machines. 

The first efficient learning procedure for large-scale Boltzmann machines 
used an extremely limited architecture, first proposed in Smolensky (1986), 
that was designed to make inference easy. A restricted Boltzmann machine 
(RBM) has a layer of visible units and a layer of hidden units with no 
connections between the hidden units. The lack of connections between 
hidden units eliminates many of the computational properties that make 
general Boltzmann machines interesting, but makes it easy to compute the 
data-dependent statistics exactly, because the hidden units are independent 
given a data vector. If connections between visible units are also prohibited. 
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the data-independent statistics can be estimated by starting Markov chains 
at hidden states that were inferred from training vectors, and alternating 
between updating all of the visible units in parallel and updating all of the 
hidden units in parallel (Hinton, 2002). It is hard to compute how many 
alternations (half-steps) of these Markov chains are needed to approach the 
stationary distribution, and it is also hard to know how close this approach 
must be for learning to make progress toward a better model. It is tempting 
to infer that if the learning works, the Markov chains used to estimate the 
data-independent statistics must be close to equilibrium, but it turns out 
that this is quite wrong. 

Empirically, learning usually works quite well if the alternating Gibbs 
sampling is run for only one full step starting from the sampled binary 
states of the hidden units inferred from a data vector (Hinton, 2002). This 
gives very biased estimates of the data-independent statistics, but it greatly 
reduces the variance in the estimated difference between data-dependent 
and data-independent statistics (Williams & Agakov, 2002), especially when 
using mini-batch learning on large data sets. Much of the sampling error in 
the data-dependent statistics caused by using a small mini-batch is elimi- 
nated because the estimate of the data-independent statistics suffers from 
a very similar sampling error. The reduced variance allows a much higher 
learning rate. Instead of viewing this learning procedure as a gross approx- 
imation to maximum likelihood learning, it can be viewed as a much better 
approximation to minimizing the difference of two divergences (Hinton, 
2002), and so it is called contrastive divergence (CD) learning. The quality 
of the learned model can be improved by using more full steps of alternat- 
ing Gibbs sampling as the weights increase from their small initial values 
(Carreira-Perpignan & Hinton, 2005), and with this modification, CD learn- 
ing allows RBMs with millions of parameters to achieve state-of-the art 
performance on a large collaborative filtering task (Salakhutdinov, Mnih, & 
Hinton, 2007). 1 

The architectural limitations of RBMs can be overcome by using them 
as simple learning modules that are stacked to form a deep, multilayer 
network. After training each RBM, the activities of its hidden units, when 
they are being driven by data, are treated as training data for the next RBM 
(Hinton, Osindero, & Teh, 2006; Hinton & Salakhutdinov, 2006). However, if 
multiple layers are learned in this greedy, layer-by-layer way, the resulting 
composite model is not a multilayer Boltzmann machine (Hinton et al., 
2006). It is a hybrid generative model called a deep belief net that has 
undirected connections between its top two layers and downward-directed 
connections between all adjacent lower layers. 


1 The performance is comparable with the best other single models, such as probabilis- 
tic matrix factorization. By averaging many models, it is possible to do better, and the 
two systems with the best performance on Netflix use multiple RBMs among the many 
models that are averaged. 
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In this article, we present a fairly efficient learning procedure for fully 
general Boltzmann machines. To estimate the data-dependent statistics, we 
use mean-field variational inference and rely on the learning to make the 
true posterior distributions be close to the factorial distributions assumed 
by the mean-field approximation. To estimate the data-independent statis- 
tics, we use a relatively small number of persistent Markov chains and rely 
on a subtle interaction between the learning and the Markov chains to al- 
low a small number of slow mixing chains to sample quickly from a highly 
multimodal energy landscape. For both sets of statistics, the fact that the pa- 
rameters are changing is essential for making the estimation methods work. 

We then show how to make our learning procedure for general Boltz- 
mann machines considerably more efficient for deep Boltzmann machines 
(DBMs) that have many hidden layers but no connections within each layer 
and no connections between nonadjacent layers. The weights of a DBM 
can be initialized by training a stack of RBMs, but with a modification 
that ensures that the resulting composite model is a Boltzmann machine 
rather than a deep belief net (DBN). This pretraining method has the added 
advantage that it provides a fast, bottom-up inference procedure for ini- 
tializing the mean-field inference. We use the MNIST and NORB data sets 
to demonstrate that DBMs learn very good generative models of images 
of handwritten digits and 3D objects. Although this article is primarily 
about learning generative models, we also show that the weights learned 
by these models can be used to initialize deep feedforward neural networks. 
These feedforward networks can then be fine-tuned using backpropagation 
to give much better discriminative performance than randomly initialized 
networks. 

2 Boltzmann Machines 

A Boltzmann machine (BM) is a network of symmetrically coupled stochas- 
tic binary units. It contains a set of visible units v e {0. 1} V and a set of 
hidden units h e {0, 1 \ lA (see Figure 1, left panel) that learn to model higher- 
order correlations between the visible units. The energy of the state {v, h} is 
defined as 


E(v, h; 6) = — v T Wh - ^v T Lv - ^h T /h, (2.1) 

where 6 = {W, L, /} are the model parameters. 2 W, L, / represent visible- 
to-hidden, visible-to-visible, and hidden-to-hidden symmetric interaction 


2 We have omitted the bias terms for clarity of presentation. Biases are equivalent to 
weights on a connection to a unit whose state is fixed at 1, so the equations for their 
derivatives can be inferred from the equations for the derivatives with respect to weights 
by simply setting the state of one of the two units to 1. 
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General Boltzmann Machine 

J 



Restricted Boltzmann Machine 



Figure 1 : (Left) A general Boltzmann machine. The top layer represents a vector 
of stochastic binary hidden variables, and the bottom layer represents a vector 
of stochastic binary visible variables. (Right) A restricted Boltzmann machine 
with no hidden-to-hidden or visible-to-visible connections. 

terms, respectively. The diagonal elements of L and / are set to 0. The 
probability that the model assigns to a visible vector v is 

P(v;6>) = P {V ' 6) = — — Vexp(-E(v, h: 6>)>, (2.2) 

2(9) 2(9) ^ F 

h 

Z(6) = '^2^2exp(-E(v,h;6)), (2.3) 

v h 


where P* denotes unnormalized probability and 2(6) is the partition func- 
tion. The conditional distributions over hidden and visible units are given 
by 

p(hj = 1 1 v, h _j)=g j J2 + Y^IjmK ) - ( 2 - 4 ) 

V ' m ^i ) 

p(v, = l|h, V_ f ) =g ^ Wfjhj + L tk v^j , (2.5) 

where g(x) = 1/(1 + exp(— x)) is the logistic function and x_- denotes a 
vector x but with x, : omitted. The parameter updates, originally derived by 
Hinton and Sejnowski (1983), that are needed to perform gradient ascent in 
the log likelihood can be obtained from equation 2.2, 

AW = a(lEp [vh T ] — E p [vh T ]), 

AL = a('Ep [vv T ] — E p [vv T ]V 

\ data model / 

AJ = a(E P [hh T ]-E P [hh T ]V 

\ data model / 
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where a is a learning rate. E p [•], the data-dependent term, is an ex- 
pectation with respect to the completed data distribution P data (h, v; 9) = 
P(h|v; 9)P ia t a (v), with P data (v) = ^ £„ <5( v ~ v ") representing the empiri- 
cal distribution, and E p [•], the data-independent term, is an expectation 

model 

with respect to the distribution defined by the model (see equation 2.2). 

Exact maximum likelihood learning in this model is intractable. The 
exact computation of the data-dependent expectation takes time that is 
exponential in the number of hidden units, whereas the exact computation 
of the model's expectation takes time that is exponential in the number of 
hidden and visible units. 

Setting both / = 0 and L = 0 recovers the restricted Boltzmann machine 
(RBM) model (see Figure 1, right panel). Setting only the hidden-to-hidden 
connections / = 0 recovers a semirestricted Boltzmann machine (Osindero 
& Hinton, 2008) in which inferring the states of the hidden units given the 
visible states is still very easy but learning is more complicated because 
it is no longer feasible to infer the states of the visible units exactly when 
reconstructing the data from the hidden states. 

2.1 A Stochastic Approximation Procedure for Estimating the Data- 
Independent Statistics. Markov chain Monte Carlo (MCMC) methods be- 
longing to the general class of stochastic approximation algorithms of the 
Robbins-Monro type (Younes, 1989; Robbins & Monro, 1951) can be used 
to approximate the data-independent statistics (Younes, 1999; Neal, 1992; 
Yuille, 2004; Tieleman, 2008). To be more precise, let us consider the follow- 
ing canonical form of the exponential family associated with the sufficient 
statistics vector <t>: 

P(X;0) = ^P) eXp(0T<t>(X)) - (Z6) 

Given a set of N independent and identically distributed (i.i.d.) training 
examples X = {x 1 , . . . , x N }, the derivative of the log likelihood with respect 
to parameter vector 9 takes the form 


3 logP(X; 9) 
3 9 


1 

N 


N 

V<b(x")-E p [<D(x)]. 

* J 1 model 

n = 1 


(2.7) 


The idea behind learning parameter vector 9 using stochastic approxima- 
tion is straightforward. Let 9 l and x f be the current parameters and the state. 
Then x f and 9 ' are updated sequentially as follows: 

• Given x', a new state x ,+1 is sampled from a transition operator 
T g , (x f+1 x f ) that leaves P(-; 9 l ) invariant (e.g., a Gibbs sampler). 
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Algorithm 1: Stochastic Approximation Algorithm for a Fully Visible Boltzmann 

Machine. 

1. Given a data set A = {x 1 ,...,x lv }. Randomly initialize 0° and M sample particles 

2. for / .= (): 7' (number of iterations) do 

3. for i = 1 : M (number of parallel Markov chains) do 

4. Sample x t+1 ’* given x 4 ’ 1 using transition operator Tgt (x i+1,! <— x t,! ). 

5. end for 

6. Update: 6 t+1 = 9* + a t Eli *(*") - M Eli <J>(* t+1 ' m )] . 

7. Decrease qj. 

8 . end for 


• A new parameter 9 t+1 is then obtained by replacing the intractable 
data -independent statistics E p [ (x ) ] by a point estimate (x f +1 ) . 

model 

In practice, we typically maintain a set of M sample points {x f l , . . . , x i,M }, 
which we often refer to as sample particles. In this case, the in- 
tractable data-independent statistics are replaced by the sample averages 
Vm Y2m = l < t > (x f+1, " z ). The procedure is summarized in algorithm 1. 

The standard proof of convergence of these algorithms relies on the fol- 
lowing basic decomposition. First, the gradient of the log-likelihood func- 
tion takes the form 


S(9) 


3 logP(X; 9) 
3 9 


1 

N 


N 

V0(x")-E p [<t>(x)]. 

‘ J 1 model 

n = 1 


( 2 . 8 ) 


The parameter update rule then takes the following form: 


9 l+1 =9 t +u t 


1 N 1 M 


n= 1 


711=1 


= 9 l + a t S(9 t ) +a t 


M 


% [4»(x)]- - V<h(x f+1 ’ m ) 

1 model ]\/l ‘ J 


— 9 ^ + q^S((^) + 


(2.9) 


The first term is the discretization of the ordinary differential equation 
9 = S(9). The algorithm is therefore a perturbation of this discretization 
with the noise term e t+1 . The proof then proceeds by showing that the noise 
term is not too large. 
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Precise sufficient conditions that ensure almost sure convergence to an 
asymptotically stable point of 9 = S(9) are given in Younes (1989, 1999) and 
Yuille (2004). One necessary condition requires the learning rate to decrease 
with time, so that a t = o o and a t < This condition can, for 
example, be satisfied simply by setting a t = a/(b + t) for positive constants 
a > 0, b > 0. Other conditions ensure that the speed of convergence of the 
Markov chain, governed by the transition operator T g/ does not decrease 
too fast as 0 tends to infinity, and that the noise term e, in the update of 
equation 2.9 is bounded. Typically, in practice, the sequence |P f | is bounded, 
and the Markov chain, governed by the transition kernel T e , is ergodic. 
Together with the condition on the learning rate, this ensures almost sure 
convergence of the stochastic approximation algorithm to an asymptotically 
stable point of 6 = S(9). 

Informally the intuition behind why this procedure works is the fol- 
lowing. As the learning rate becomes sufficiently small compared with the 
mixing rate of the Markov chain, this persistent chain will always stay 
very close to the stationary distribution, even if it is run for only a few 
MCMC steps per parameter update. Samples from the persistent chain will 
be highly correlated for successive parameter updates, but if the learning 
rate is sufficiently small, the chain will mix before the parameters have 
changed enough to significantly alter the value of the estimator. 

The success of learning relatively small Boltzmann machines (Neal, 1992) 
seemed to imply that the learning rate was sufficiently small to allow the 
chains to stay close to equilibrium as the parameters changed. Recently, 
however, this explanation has been called into question. After learning an 
RBM using persistent Markov chains for the data-independent statistics, we 
tried sampling from the RBM and discovered that even though the learning 
had produced a good model, the chains mixed extremely slowly. In fact, 
they mixed so slowly that the appropriate final learning rate, according 
to the explanation above, would have been smaller by several orders of 
magnitude than the rate we actually used. So why did the learning work? 

Tieleman and Hinton (2009) argue that the fact that the parameters are 
being updated using the data-independent statistics gathered from the per- 
sistent chains means that the mixing rate of the chains with their parameters 
fixed is not what limits the maximum acceptable learning rate. Consider, 
for example, a persistent chain that is stuck in a deep local minimum of 
the energy surface. Assuming that this local minimum has very low prob- 
ability under the posterior distributions that are used to estimate the data- 
dependent statistics, the effect of the learning will be to raise the energy of 
the local minimum. After a number of weight updates, the persistent chain 
will escape from the local minimum not because the chain has had time 
to mix but because the energy landscape has changed to make the local 
minimum much less deep. The learning causes the persistent chains to be 
repelled from whatever state they are currently in, and this can cause slow- 
mixing chains to move to other parts of the dynamic energy landscape much 
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faster than would be predicted by the mixing rate with static parameters. 
Welling (2009) has independently reported a closely related phenomenon 
that he calls "herding." 

Recently Tieleman (2008), Salakhutdinov and Hinton (2009a), Salakhut- 
dinov (2009), and Desjardins, Courville, Bengio, Vincent, and Delalleau 
(2010) have shown that this stochastic approximation algorithm, also 
termed persistent contrastive divergence, performs well compared to con- 
trastive divergence at learning good generative models in RBMs. Although 
the allowable learning rate is much higher than would be predicted from 
the mixing rate of the persistent Markov chains, it is still considerably lower 
than the rate used for contrastive divergence learning because the gradient 
estimate it provides has lower bias but much higher variance. The vari- 
ance is especially high when using online or mini-batch learning rather 
than full batch learning. With contrastive divergence learning, the error 
in the data-dependent statistics introduced by the fact that a mini-batch 
has sampling error is approximately cancelled by approximating the data- 
independent statistics using short Markov chains that start at the data in 
the mini-batch and do not have time to move far away. With persistent 
contrastive divergence, this strong correlation between the sampling errors 
in the data-dependent and data-independent statistics is lost. 


2.2 A Variational Approach to Estimating the Data-Dependent Statis- 
tics. Persistent Markov chains are less appropriate for estimating the data- 
dependent statistics, especially with mini-batch learning on large data sets. 
Fortunately, variational approximations work well for estimating the data- 
dependent statistics. Given the data, it is typically quite reasonable for the 
posterior distribution over latent variables to be unimodal, especially for 
applications like speech and vision where normal data vectors really do 
have a single correct explanation and the data are rich enough to allow that 
explanation to be inferred using a good generative model. 

In variational learning (Zemel, 1993; Hinton & Zemel, 1994; Neal & 
Hinton, 1998; Jordan et al., 1999), the true posterior distribution over latent 
variables P(h|v; 9) for each training vector v is replaced by an approxi- 
mate posterior Q(h|v; /z), and the parameters are updated to maximize the 
variational lower bound on the log likelihood. 


logP(v; 9) > ^ Q(h|v; /z) logP(v, h; 6) + H(Q) 

h 

>logP(v; 9) - KL[Q(h|v; /z)||P(h|v; 6>)] , (2.10) 


where H(-) is the entropy functional. 

Variational learning has the nice property that in addition to trying to 
maximize the log likelihood of the training data, it tries to find parameters 
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that minimize the Kullback-Leibler divergence between the approximating 
and true posteriors. Making the true posterior approximately unimodal, 
even if it means sacrificing some log likelihood, could be advantageous 
for a system that will use the posterior to control its actions. Having mul- 
tiple alternative representations of the same sensory input increases the 
likelihood compared with a single explanation of the same quality, but it 
makes it more difficult to associate an appropriate action with that sensory 
input. Variational inference that uses a factorial distribution to approxi- 
mate the posterior helps to eliminate this problem. During learning, if the 
posterior given a training input vector is multimodal, the variational in- 
ference will lock onto one mode, and learning will make that mode more 
probable. Our learning algorithm will therefore tend to find regions in 
the parameter space in which the true posterior is dominated by a single 
mode. 

For simplicity and speed, we approximate the true posterior using a 
fully factorized distribution (i.e., the naive mean-field approximation), 
Q(h; fi) = ri/li where cj{hj = 1) = /x • and U is the number of hidden 
units. The lower bound on the log probability of the data takes the following 
form: 


The learning proceeds by first maximizing this lower bound with respect 
to the variational parameters // for fixed 9, which results in the mean-field 
fixed-point equations: 


This is followed by applying stochastic approximation to update model 
parameters 9, which is summarized in algorithm 2. 

We emphasize that variational approximations should not be used for es- 
timating the data-independent statistics in the Boltzmann machine learning 
rule, as attempted in Galland (1991), for two separate reasons. First, a facto- 
rial approximation cannot model the highly multimodal, data-independent 
distribution that is typically required. Second, the minus sign causes the 
parameters to be adjusted so that the true model distribution becomes as 
different as possible from the variational approximation. 


logP(v; 9) > * 



^[(Uylog/Uy + (1 - flj) log (1 - IAj)\. 



( 2 . 11 ) 
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Algorithm 2: Learning Procedure for a General Boltzmann Machine. 

1. Given: a training set of N binary data vectors {v}^_ 1 , and M, the number of persistent 
Markov chains (i.e., particles). 

2. Randomly initialize parameter vector 9 0 and M samples: {v 0,1 , h 0,1 }, {v 0,M , h 0,M }. 

3. for t = 0 to T (number of iterations) do 

4. // Variational Inference: 

5. for each training example \ n ,n = 1 to N do 

6. Randomly initialize /i and run mean-field updates until convergence: 


7. Set fi n — fj,. 

8. end for 

9. // Stochastic Approximation: 

10. for each sample m = 1 to M (number of persistent Markov chains) do 

11. Sample (v £+1,m , h £+1,m ) given (v £,m , h i,m ) by running a Gibbs sampler 
(see equations 2.4, 2.5) 

12. end for 

13. // Parameter Update : 


17. Decrease at. 

18. end for 


3 Learning Deep Boltzmann Machines 

Algorithm 2 can learn Boltzmann machines with any pattern of connectivity 
between the emits, but it can be made particularly efficient in deep Boltz- 
mann machines that have multiple hidden layers but have connections only 
between adjacent layers, as shown in Figure 2. Deep Boltzmann machines 
are interesting for several reasons. First, like deep belief networks, DBMs 
have the ability to learn internal representations that capture very complex 
statistical structure in the higher layers. As has already been demonstrated 
for DBNs, this is a promising way of solving object and speech recognition 
problems (Bengio, 2009; Bengio & LeCun, 2007; Flinton et al., 2006; Dahl, 
Ranzato, Mohamed, & Hinton, 2010; Mohamed, Dahl, & Hinton, 2012). 
High-level representations can be built from a large supply of unlabeled 
data, and a much smaller supply of labeled data can then be used to fine- 
tune the model for a specific discrimination task. Second, again like DBNs, 
if DBMs are learned in the right way, there is a very fast way to initialize 
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Deep Belief Network Deep Boltzmann Machine 




Figure 2: (Left) Deep belief network (DBN). The top two layers form an 
undirected graph, and the remaining layers form a belief net with directed, 
top-down connections (Right) Deep Boltzmann machine (DBM), with both 
visible-to-hidden and hidden-to-hidden connections but no within-layer conn- 
ections. All the connections in a DBM are undirected. 


the states of the units in all layers by simply doing a single bottom-up pass 
using twice the weights to compensate for the initial lack of top-down feed- 
back. Third, unlike DBNs and many other models with deep architectures 
(Ranzato, Huang, Boureau, & LeCun, 2007; Vincent, Larochelle, Bengio, 
& Manzagol, 2008; Serre, Oliva, & Poggio, 2007), the approximate in- 
ference procedure, after the initial bottom-up pass, can incorporate top- 
down feedback. This allows DBMs to use higher-level knowledge to 
resolve uncertainty about intermediate-level features, thus creating bet- 
ter data-dependent representations and better data-dependent statistics for 
learning. 3 

Let us consider a three-hidden-layer DBM, as shown in Figure 2 (right), 
with no within-layer connections. The energy of the state {v, h (1) , h <2) , h (31 } 
is defined as 


E(v, h (1) , h (2) , h (3) ; 6 ) = -v T W (1) h (1) - h (llT W (2) h (2) - h (2)T W (3) h <3) , 

(3.1) 


3 For many learning procedures, there is a trade-off between the time taken to infer 
the states of the latent variables and the number of weight updates required to learn a 
good model. For example, an autoencoder that uses noniterative inference requires more 
weight updates than an autoencoder that uses iterative inference to perform a look-ahead 
search for a code that is better at reconstructing the data and satisfying penalty terms 
(Ranzato, 2009). 
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where 9 = {W (1) , W (2) , W (3) } are the model parameters, representing 
visible-to-hidden and hidden-to-hidden symmetric interaction terms. 

The probability that the model assigns to a visible vector v is 

P(v;6>) = — — exp (— E(v, h (1) , h (2) , h (3) ; 9)). (3.2) 

hW.hM.h® 


The conditional distributions over the visible and the three sets of hidden 
units are given by logistic functions: 


paf = nv, h< 2 >) = g J2 wVvi + J2 w]X 

\ i m 

X = l|h (1) , h <3) ) =g (£ + J2 W M h i 3) 

P (h^ = X)=§(j2 w LX 


p(«i = l|h (1) )=g 


(3.3) 

(3.4) 

(3.5) 

(3.6) 


The learning procedure for general Boltzmann machines described above 
can be applied to DBMs that start with randomly initialized weights, but 
it works much better if the weights are initialized sensibly. With small 
random weights, hidden units in layers that are far from the data are very 
underconstrained, so there is no consistent learning signal for their weights. 
With larger random weights, the initialization imposes a strong random 
bias on the feature detectors learned in the hidden layers. Even when the 
ultimate goal is some unknown discrimination task, it is much better to bias 
these feature detectors toward ones that form a good generative model of 
the data. We now describe how this can be done. 


3.1 Greedy Layerwise Pretraining of DBNs. Hinton et al. (2006) intro- 
duced a greedy, layer-by-layer unsupervised learning algorithm that con- 
sists of learning a stack of RBMs one layer at a time. After greedy learning, 
the whole stack can be viewed as a single probabilistic model called a deep 
belief network. Surprisingly, this composite model is not a deep Boltzmann 
machine. The top two layers form a restricted Boltzmann machine, but the 
lower layers form a directed sigmoid belief network (see Figure 2, left). 
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After learning the first RBM in the stack, the generative model can be 
written as 


P(v; 9) = ^P(h (1) ; W (1) )P(v|h a) ; W a) ), (3.7) 

h<D 

where P(h (1) ; W (1) ) = ^ v P(h (1) , v; W (1) ) is a prior over h (1) that is implic- 
itly defined by W (1) . Using the same parameters to define both the prior 
over h (1) and the likelihood term P(v|h (11 ) seems like an odd thing to do for 
those who are more familiar with directed graphical models, but it makes 
inference much easier and it is only a temporary crutch: the prior over h (l) 
defined by W (1) will be thrown away and replaced by a better prior defined 
by the weights, W (2) , of the next RBM in the stack. 

The second RBM in the stack attempts to learn a better overall model by 
leaving P(v|h (1) ; W (1) ) fixed and replacing P(h (1) ; W (1) ) by 


P(h (1) ; W (2) ) = ^]P(h (1) ,h (2) ; W (2) ), 

h< 2 > 

where W <2) is initialized at W (1)T and then improved by following the 
gradient of a variational lower bound on the log probability of the training 
data with respect to IV a> . The variational bound was first derived using 
coding arguments in Hinton and Zemel (1994). For a data set containing N 
training examples, it has the form 


N 

£log P(v"; 0)>^E Q(ha)|vll) [logP(v"|h a >; W (1) )] 

n = 1 n 


-^KL(Q(h a) |v")||P(h a) ; W (2) )) 

n 


^£ 


Y, Q(h a V)[logP(v|h«; W* 1 ’)] + H(Q) 


n Lh<D J 

+ £ £ Q(h a) |v")logP(h (1) ; W (2) ), (3.8) 

n h« 


where H{-) is the entropy functional and Q(h (1) |v) is any approximation 
to the posterior distribution over vectors in hidden layer 1 for the DBN 
containing hidden layers h (1) and h ,2) . The approximation we use is the 
true posterior over h (1) for the first RBM, Q(h (11 |v) = P(h (1) |v; W (11 ). If the 
second RBM is initialized to be the same as the first RBM but with its visible 
and hidden units interchanged, W (2) = W (1)T , then Q(h (1) |v) defines the 
DBN's true posterior over h (1) and the variational bound is tight. As soon 
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as W (2) ceases to be identical to W (1)T , Q is no longer the true posterior for 
the DBN. 

Changing W (2) affects only the last sum in equation 3.8, so maximiz- 
ing the bound, summed over all the training cases, with regard to W (2 ^ 
amounts to learning a better model of the mixture over all N training cases: 
1 Q(h (1) |v"), which we call the aggregated posterior. Each individual 

posterior of the first RBM Qfh' 1 - 1 |v" ) is factorial, but the aggregated poste- 
rior is typically very far from factorial. Changing W (2) so that the second 
RBM becomes a better model of the aggregated posterior over h (1) is then 
guaranteed to improve the variational bound for the whole DBN on the log 
likelihood of the training data. 

This argument can be applied recursively to learn as many layers of fea- 
tures as desired. Each RBM in the stack performs exact inference while 
it is being learned, but once its implicit prior over its hidden vectors 
has been replaced by a better prior defined by the higher-level RBM, the 
simple inference procedure ceases to be exact. As the stack gets deeper, the 
simple inference procedure used for the earlier layers can be expected to 
become progressively less correct. Nevertheless, each time a new layer is 
added and learned, the variational bound for the deeper system is better 
than the bound for its predecessor. When a third hidden layer is added, for 
example, the bound of equation 3.8 is replaced by a bound in which the last 
sum, 


E E Q(h (1) |v") logP(h (1) ; W (2) ), (3.9) 

n hM 

is replaced by 

E E Q (h (1) I v" ) (Eq ( i,<2) | h d) ) [log P(h (1 ’ | h (2) ; W ,2) )] 

n hW 

-KL(Q(h (2, |h (1) )l|P(h <2) ; W (3) ))). (3.10) 

When the second RBM is learned, the log probability of the training data also 
improves because the bound starts off tight. For deeper layers, we are still 
guaranteed to improve the variational lower bound on the log probability 
of the training data, though the log probability itself can decrease. For 
these deeper layers, the bound does not start off tight and could therefore 
become tighter, allowing the log probability to decrease even though the 
bound increases. 

The improvement of the bound is guaranteed only if each RBM in the 
stack starts with the same weights as the previous RBM and follows the gra- 
dient of the log likelihood, using the posterior distributions over the hidden 
units of the previous RBM as its data. In practice, we violate this condition by 
using gross approximations to the gradient such as contrastive divergence. 
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The real value of deriving the variational bound is to allow us to under- 
stand why it makes sense to use the aggregated posterior distributions of 
one RBM as the training data for the next RBM and why the combined 
model is a deep belief net rather than a deep Boltzmann machine. 

3.2 Greedy Layerwise Pretraining of DBMs. Although the simple way 
of stacking RBMs leads to a deep belief net, it is possible to modify the pro- 
cedure so that stacking produces a deep Boltzmann machine. We start by 
giving a crude intuitive argument about how to combine RBMs to get a 
deep Boltzmann machine. We then show that a modification of the intu- 
itively derived method for adding the top layer is guaranteed to improve a 
variational bound. We also show that the procedure we use in practice for 
adding intermediate layers fails to achieve the property that is required for 
guaranteeing an improvement in the variational bound. 

After training the second-layer RBM in a deep belief net, there are 
two ways of computing a factorial approximation to the true poste- 
rior P(h (i:> |v; Wfl'Wffl). The obvious way is to ignore the second-layer 
RBM and use the P(h ( 1 ) |v; W (1> ) defined by the first RBM. An alternative 
method is to first sample h (1) from P(h ( 1 ) |v; W (1> ), then sample h (2) from 
P(h ( 2 ) |h (1) ; W (2) ), and then use the P(h ( 1 ) |h (2) ; W (2) ) defined by the second 
RBM . 4 The second method will tend to overemphasize the prior for h (1) 
defined by W <2> , whereas the first method will tend to underemphasize this 
prior in favor of the earlier prior defined by W (1) that it replaced. 

Given these two different approximations to the posterior, it would be 
possible to take a geometric average of the two distributions. This can be 
done by first performing a bottom-up pass to infer h (2) then using V 2 W® 
and V 2 W (2) to infer h (1) from both v and h (2) . Notice that h (2) is inferred 
from v, so it is not legitimate to sum the full top-down and bottom-up 
influences. This would amount to double-counting the evidence provided 
by v and would give a distribution that was much too sharp. Experiments 
with trained DBNs confirm that averaging the top-down and bottom-up 
inputs works well for inference and adding them works badly. 

This reasoning can be extended to a much deeper stack of greedily trained 
RBMs. The initial bottom-up inference that is performed in a DBN can be 
followed by a stage in which all of the weights are halved and the states 
of the units in the intermediate layers are resampled by summing the top- 
down and bottom-up inputs to a layer. If we alternate between resampling 
the odd-numbered layers and resampling the even-numbered layers, this 
corresponds to alternating Gibbs sampling in a deep Boltzmann machine 
with the visible units clamped. So after learning a stack of RBMs, we can 
either compose them to form a DBN or halve all the weights and compose 


4 The sampling noise in the second method can be reduced by using a further approx- 
imation in which the sampled binary values are replaced by their probabilities. 
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Figure 3: Pretraining a DBM with three hidden layers consists of learning a 
stack of RBMs that are then composed to create a DBM. The first and last RBMs 
in the stack need to be modified by using weights that are twice as big in one 
direction. 


them to form a DBM. Moreover, given the way the DBM was created, there 
is a very fast way to initialize all of the hidden layers when given a data 
vector: simply perform a bottom-up pass using twice the weights of the 
DBM to compensate for the lack of top-down input. 

There is an annoying problem with this method of pretraining a DBM. 
For the intermediate layers, using weights that are half of the weights 
learned by the individual RBMs seems fine because when the RBMs are 
combined, it can be viewed as taking the geometric mean of the bottom-up 
and top-down models, but for the visible layer and the top layer, it is not 
legitimate because they receive input only from one other layer. Both the 
top and the bottom layers need to be updated when estimating the data- 
independent statistics for the DBM, and we cannot use weights that are 
bigger in one direction than the other because this does not correspond to 
Gibbs sampling in any energy function. So we need to use a special trick 
when pretraining the first and last RBMs in the stack. 

For the first RBM, we constrain the bottom-up weights to be twice the top- 
down weights during the pretraining, as shown in Figure 3. This means that 
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we can halve the bottom-up weights without halving the top-down weights 
and still be left with symmetric weights. Conversely, for the last RBM in 
the stack, we constrain the top-down weights to be twice the bottom-up 
weights. Now we have a stack of RBMs that we can convert to a DBM by 
halving all but the first layer top-down weights and the last layer bottom-up 
weights. 5 

Greedily pretraining the weights of a DBM in this way serves two pur- 
poses. First, it initializes the weights to reasonable values. Second, it ensures 
that there is a fast way of performing approximate inference by a single 
bottom-up pass using twice the weights for all but the top-most layer. This 
eliminates the need to store the hidden states that were inferred last time 
a training case was used (Neal, 1992) or to use simulated annealing from 
random initial states of the hidden units (Hinton & Sejnowski, 1983). This 
fast approximate inference is used to initialize the mean-field, iterative in- 
ference, which then converges much faster than mean field with random 
initialization. Since the mean-field inference uses real-valued probabilities 
rather than sampled binary states, we also use probabilities rather than 
sampled binary states for the initial bottom-up inference. 

3.3 A Variational Bound for Greedy Layerwise Pretraining of a DBM 
with Two Hidden Layers. The explanation of DBM pretraining given in 
the previous section is motivated by the need to end up with a deep network 
that has symmetric weights between all adjacent pairs of layers. However, 
unlike the pretraining of a DBN, it lacks proof that each time a layer is 
added to the DBM, the variational bound for the deeper DBM is better than 
the bound for the previous DBM. We now give an explanation of why the 
method for training the first RBM in the stack works, and we show that with 
a slight modification, the proposed method for training the second layer 
RBM is a correct way of improving a variational bound. The main point 
of this exercise is to gain a better understanding of how the pretraining 
procedures for DBNs and DBMs are related. The apparently unproblematic 
method for pretraining the intermediate layers using a symmetric RBM is 
actually more problematic than it seems, and we discuss this in section 3.4. 

The basic idea for pretraining a DBM is to start by learning a model 
in which the prior over hidden vectors, p(h (1) ; W (1) ), is the normalized 
product of two identical distributions. Then one of these distributions is 
discarded and replaced by the square root of a better prior p(h (1) ; W <2) ) that 
has been trained to fit a good approximation to the aggregated posterior of 
the first model. 


5 Of course, when we constrain the weights of an RBM to have different magnitudes 
in the two directions, the usual rules for updating the states of the units no longer corre- 
spond to alternating Gibbs sampling in any energy function, but the one-step contrastive 
divergence learning still works well, for reasons that will be explained later. 
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Figure 4: Pretraining a deep Boltzmann machine with two hidden layers, 
(a) The DBM with tied weights is trained to model the data using one-step 
contrastive divergence (see Figure 5). (b) The second hidden layer is removed, 
(c) The second hidden layer is replaced by part of the RBM that has been trained 
on the aggregated posterior distribution at the first hidden layer of the DBM in 
panel a. (d) The resulting DBM with a modified second hidden layer. 


Figure 5 : Pretraining the DBM with two hidden layers shown in Figure 4a using 
one-step contrastive divergence. The second hidden layer is initialized to be the 
same as the observed data. The units in the first hidden layer have stochastic 
binary states, but the reconstructions of both the visible and second hidden 
layer use the unsampled probabilities, so both reconstructions are identical. The 
pairwise statistics for the visible layer and the first hidden layer are therefore 
identical to the pairwise statistics for the second hidden layer and the first 
hidden layer, and this is true for both the data and the reconstructions. 

Consider the simple DBM in Figure 4a that has two hidden layers and 
tied weights. If we knew what initial state vector to use for h (2) , we could 
train this DBM using one-step contrastive divergence with mean-field re- 
constructions of both the visible states and the states of the top layer, as 
shown in Figure 5. So we simply set the initial state vector of the top layer 
to be equal to the data, v. Provided we use mean-field reconstructions for 
the visible units and the top-layer units, one-step contrastive divergence is 
then exactly equivalent to training an RBM with only one hidden layer but 
with bottom-up weights that are twice the top-down weights, as prescribed 


h( 2 ) = v 



v 
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Algorithm 3: Greedy Pretraining Algorithm for a Deep Boltzmann Machine. 

1. Train the first-layer “RBM” using one-step contrastive divergence learning with mean-held 
reconstructions of the visible vectors. During the learning, constrain the bottom-up weights, 
2 IpW, to be twice the top-down weights, Ifd 1 ). 

2. Freeze 2 IP 9) that defines the first layer of features and use samples h*P from 
P(hW|v; 2 fP( 1 )) as the data for training the second RBM. This is a proper RBM with 
weights 2 that are of the same magnitude in both directions. It is also trained us- 
ing one-step contrastive divergence learning with mean-held reconstructions of its visible 
vectors. 

3. Freeze 2W^ that dehnes the second layer of features and use the samples lb 2 ) 
from P(h( 2 )|v; 21PP1. 2ip( 2 ^) as the data for training the next RBM in the same way 
as the previous one. 

4. Proceed recursively up to layer L — 1. 

5. Train the top-level RBM using one-step contrastive divergence learning with mean-held 
reconstructions of its visible vectors. During the learning, constrain the bottom-up weights, 
W( L \ to be half the top-down weights, 2W^ L \ 

6. Use the weights {W^\ W^ 2 \ W^\ ..., W^} to compose a deep Boltzmann machine. 


in algorithm 3. This way of training the simple DBM with tied weights is un- 
likely to maximize the likelihood of the weights, but in practice it produces 
surprisingly good models that reconstruct the training data well. 

After pretraining, we can use 2 W (1) to compute an approximate posterior 
Q(h (1) |v) for h (1) from v, and using this approximate posterior, we can get 
a variational lower bound on the log probability (see equation 3.8), that the 
simple DBM in Figure 4a assigns to the training data: 

lQ g p ( v n } > E Q(hW |v„) Dog p I h (1) ; w (1) )] 

n n 

-^KL(Q(h (1 »|v n )||P(h (1) ;W (1) )). (3.11) 


We now show how to improve this variational bound. 

The model's marginal distribution over h (1) is the product of two identi- 
cal distributions — one defined by an RBM composed of h (1) and v and the 
other defined by an identical RBM composed of h (1) and h (2) : 

P(h U> ; W CD) = _L_ ^ e vwJ ^ e h <2 > T W<'>h<^ , (3.1.2) 
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where Z(\N ( '' ) ) is the normalizing constant. 6 The idea is to keep one of 
these two RBMs and replace the other by the square root of a better prior 
P(h (1) ; W <2) ). To do so, we train another RBM with two sets of hidden 
units and tied weights, as shown in Figure 4c, to be a better model of the 
aggregated variational posterior ^ Yin Q(h (l, |v n ; W (1> ) of the first model 
(see Figure 4a). If the tied weights are initialized at = W (1)T , then the 
higher-level RBM has exactly the same prior over h (1) as the original DBM. 
If the RBM is trained by following the gradient of the log likelihood, 7 any 
change in the weights will ensure that 


Y KL(Q(h (1) |v„; W (1) )||P(h (1) ; W (2) )) 


< Y KL(Q(h (1) |v„; W (1) )||P(h (1) ; W (1) )). (3.13) 

n 


Similar to equation 3.12, the distribution over h (1) defined by the second- 
layer RBM is also the product of two identical distributions, one for each set 
of hidden units. This implies that taking a square root amounts to simply 
keeping one such distribution. 

Once the two RBMs are composed to form a two-hidden-layer DBM 
model (see Figure 4d), the marginal distribution over h' 1 - 1 is the geometric 
mean of the two probability distributions, P(h fl) ; W (1) ), P(h (1) ; W (2) ), de- 
fined by the first- and second-layer RBMs (i.e., the renormalized pairwise 
products of the square roots of the two probabilities for each event): 

P(h “ ■ w "’ w<2>) = (£« l " T,w ")• 

(3.14) 

The variational lower bound of equation 3.11 improves because replacing 
half of the prior by a better model reduces the Kullback-Leibler divergence 
from the variational posterior, as shown in the appendix. 

Due to the convexity of asymmetric divergence, this is guaranteed to 
improve the variational bound of the training data by at least half as much 


6 The biases learned for h (1) are shared equally between the two RBMs. 

7 In practice, we depart from maximum likelihood training is several ways: we use 
one-step contrastive divergence to train the second-layer RBM; we use mean-field recon- 
structions of its visible units; we do not sample its two sets of hidden units independently. 
But the same issue arises with the variational bound for DBNs. That bound also requires 
proper maximum likelihood training of the higher-level RBMs. Fortunately, during pre- 
training, we typically terminate the learning of each RBM long before it converges, and 
in this early regime, contrastive divergence learning is almost always improving the 
likelihood, which is all we need. 
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as fully replacing the original prior. It is also guaranteed to loosen the 
variational bound by at most half as much as fully replacing the original 
prior, assuming that inference is still performed assuming the original prior. 

This argument shows that the apparently unprincipled hack of doubling 
the weights in one direction to cope with the "end effects" when creating a 
two-hidden-layer DBM from a stack of two RBMs can actually be justified 
as a way of improving a variational bound, except for the fact that the top 
RBM should be trained by maximum likelihood. 


3.4 Pretraining Intermediate Hidden Layers. We have not been able 
to design a way of adding intermediate hidden layers that is guaranteed to 
improve a variational bound, but the scheme we use in algorithm 3 seems 
to work well in practice. One difficulty in extending the proof is that for an 
intermediate RBM, we need to take the square root of the marginal prior 
distribution over its visible units so that we can use it to replace half of the 
previous prior of the lower RBM. We also need to take the square root of 
the marginal prior distribution over the hidden units of the intermediate 
RBM in order to allow the next RBM to replace half of this prior. 

Halving the weights of an RBM halves the energy of every joint config- 
uration of the hidden and visible units, so it takes the square root of the 
joint distribution over pairs of visible and hidden vectors, and it also takes 
the square root of the conditional distribution over visible vectors given a 
hidden vector (or vice versa), but it does not take the square root of the 
marginal prior distributions over the visible or the hidden vectors. This is 
most easily seen by considering the ratios of the probabilities of two visible 
vectors, v ( , and v ; , . Before halving the weights, their probability ratio in the 
marginal distribution over visible vectors is given by 


P(yJ _ Eh e ~ £(v “ ,hl 

P( V E h e' E(vh) 


(3.15) 


In the RBM with halved weights, all of the exponents are halved which 
takes the square root of every individual term in each sum, but this does 
not take the square root of the ratio of the sums. This argument shows 
that the apparently unproblematic idea of halving the weights of all the 
intermediate RBMs in the stack is not the right thing to do if we want to 
ensure that as each layer is added, the new variational bound is better than 
the old one. 8 


8 Removing one of the hidden groups in Figure 4c halves the expected input to the 
visible units, but it also halves the entropy of the hidden units. This halves the free energy 
of each visible vector, which takes the square root of the marginal prior distribution over 
the visible units. 
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Although we cannot guarantee that a variational bound is improved 
when adding intermediate layers, the idea that each new RBM is approxi- 
mately replacing half of the prior over the top layer of the previous model 
is a useful insight into what is happening during pretraining. It suggests, 
for example, that each new layer does less work in a DBM than it does in a 
DBN, which replaces all of the prior of the previous RBM. This in turn sug- 
gests that in the pretraining phase, DBMs are likely to get less benefit than 
DBNs from being very deep. A different method of pretraining DBMs that 
distributes the modeling work more evenly over the layers would probably 
be helpful. 

4 Evaluating Deep Boltzmann Machines 

Assessing the generalization performance of DBMs plays an important role 
in model selection, model comparison, and controlling model complexity. 
In this section we discuss two ways of evaluating the generalization capa- 
bilities of DBMs: generative and discriminative. 

4.1 Evaluating DBMs as Generative Models. We first focus on eval- 
uating generalization performance of DBMs as density models. For many 
specific tasks, such as classification or information retrieval, performance of 
DBMs can be directly evaluated (see section 4.2). More broadly, however, the 
ability of DBMs to generalize can be evaluated by computing the probability 
of held-out input vectors, which is independent of any specific application. 
An unfortunate limitation of DBMs is that the exact computation of the 
probability that the model assigns to a visible vector v is intractable. It re- 
quires both an intractable summation over the hidden vectors at all layers 
and an intractable computation of the partition function: 


Recently, Salakhutdinov and Murray (2008) showed that a Monte Carlo- 
based method, annealed importance sampling (AIS; Neal, 2001), can be 
used to efficiently estimate the partition function of an RBM. In this section, 
we show how AIS can be used to estimate the partition functions of DBMs. 
Together with variational inference, this will allow us to obtain good es- 
timates of the lower bound on the log probability of the training and test 
data. 

Suppose we have two distributions defined on some space X with prob- 
ability density functions P A (x) = P A (x) /Z A and P B (x) = P B (x) /Z B . Typically 
P A (x) is defined to be some simple distribution, with known partition func- 
tion Z A , from which we can easily draw i.i.d. samples. AIS estimates the ra- 
tio Z B j Z A by defining a sequence of intermediate probability distributions: 



hm.hcyh® 


exp (— E(v, h (1) , h (2) , h (3) ; 6)). 


(4.1) 
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P Q , P 1 , . . . , P K , with P 0 = P A and P K = P B , which satisfy P k (x) # 0 whenever 
P k+ 1 ( x > 7^ 0, k = 0, K— 1. For each intermediate distribution, we must be 
able to easily evaluate the unnormalized probability P* (x), and we must 
also be able to sample x' given x using a Markov chain transition opera- 
tor T k (x' <r- x) that leaves P k (x) invariant. One general way to define this 
sequence is to set 

P i (x)aP;(x) 1 -^P*(x)^, (4.2) 

with 0 = f) 0 < < • ■ ■ < = 1 chosen by the user. 

Using the special layer-by-layer structure of DBMs, we can derive an ef- 
ficient AIS scheme for estimating the model's partition function. Let us con- 
sider a three-hidden-layer Boltzmann Machine (see Figure 3, right) whose 
energy is defined as 

E(v, h (1) , h <2) , h (3) ; 9) = -v T W ( 1 ) h (1) - h ( 1 )T W ( 2 ) h (2) - h ( 2 )T W ( 3 ) h (3) . 

(4.3) 

By summing out the first- and the third-layer hidden units {h ( 1 ) ,h (3) }, 
we can easily evaluate an unnormalized probability P*(v, h (2) ; 9). We can 
therefore run AIS on a much smaller state space x = {v, h (2) ) with h (1) and 
h (3) analytically summed out. The sequence of intermediate distributions, 
parameterized by fi, is defined as 

P k (y, h (2) ; 9) 

= — PjTfv , h (2) ; 9) = — P k (v - h<1) > h<2) - h<3) ; e '> 

Zk Z k h a,, h (3) 

= i. Y\(l + Y\( 1 + e^«)). (4.4) 

k i i 

We gradually change (the inverse temperature) from 0 to 1, annealing 
from a simple "uniform" model P 0 to the complex deep Boltzmann machine 
model P K . 

Using equations 3.3 to 3.6, it is straightforward to derive a Gibbs transi- 
tion operator T k ([v, h (2) }' <— {v, h (2) }) that leaves P k (v, h <2) ; 9) invariant: 


= iiv, h <2) ) = g U X) K )v ‘ + E w iX 


(4.5) 


P (hM = iih« h< 3 >) = g u Y,K>? + E w ^r 


(4.6) 
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Algorithm 4: Annealed Importance Sampling Run. 


1. Given a sequence 0 = /3o < Pi < ■■■ < /3k = 1. Let x = {v, h^ 2 )}. 

2. Sample x 1 from Pq( x )- 

3. for k = 2 : K do 

4. Sample x fc given x* _1 using <— x fc_1 ) (see equations 4.5-4. 8). 

5. end for 

6. Compute importance weight using equation 4.4: 



P (/z ; < 3) = iih' 2 >)=g 


(4.7) 


m 


P(v i = ^ a) )=g\P k 'E W $ ) h? 


(4.8) 


Algorithm 4 shows a single run of AIS. Note that there is no need to compute 
the normalizing constants of any intermediate distributions. After perform- 
ing M runs of AIS, the importance weights w ® can be used to obtain an 
estimate of the ratio of partition functions: 


where Z 0 = 2 |V|+|W| is the partition function of the uniform model (e.g., 
/1 0 = 0) with |V| and \U\ denoting the total number of visible and hidden 
units, and Z K is the partition function of the DBM model (e.g., p K = 1). 

Provided K is kept large, the total amount of computation can be split 
in any way between the number of intermediate distributions K and the 
number of annealing runs M without adversely affecting the accuracy of 
the estimator. The number of AIS runs can be used to control the variance 
in the estimate of r AIS , as samples drawn from P 0 are independent: 



(4.9) 



(4.10) 


where s 2 is estimated simply from the sample variance of the importance 
weights. 
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Finally, once we obtain an estimate of the global partition function Z, 
we can estimate, for a given test case v*, the variational lower bound of 
equation 2.10: 

logP(v*; d)>-J2 Q( h l v *; aOE(v*, h; 9) + H(Q ) - log Z(0) 

h 

« - Q (h l v * ; h : + W(Q) - logi, 

h 

where we defined h = {h (1 \ h' :2) . h (3) ). For each test vector under consid- 
eration, this lower bound is maximized with respect to the variational pa- 
rameters ii using the mean-field update equations. 

Furthermore, by explicitly summing out the states of the hidden units 
{h (2) , h (3) }, we can obtain a tighter variational lower bound on the log 
probability of the test data. Of course, we can also adopt AIS to estimate 
P*(y) = EhAh«,h<3. £ *( v ’ h (1 \ h (2 \ h (3) ), and together with an estimate of 
the global partition function, we can estimate the true log probability of 
the test data. This, however, would be computationally very expensive, 
since we would need to perform a separate AIS run for each test case. As 
an alternative, we could adopt a variation of the Chib-style estimator, pro- 
posed by Murray and Salakhutdinov (2009). In the case of deep Boltzmann 
machines, where the posterior over the hidden units tends to be unimodal, 
their proposed Chib-style estimator can provide good estimates of log P* (v) 
in a reasonable amount of computer time. 

In general, when learning a deep Boltzmann machine with more than 
two hidden layers and no within-layer connections, we can explicitly sum 
out either odd or even layers. This will result in a better estimate of the 
model's partition function and tighter lower bounds on the log probability 
of the test data. 

4.2 Evaluating DBMs as Discriminative Models. After learning, the 
stochastic activities of the binary features in each layer can be replaced by de- 
terministic, real-valued probabilities, and a deep Boltzmann machine with 
two hidden layers can be used to initialize a multilayer neural network in 
the following way. For each input vector v, the mean-field inference is used 
to obtain an approximate posterior distribution Q(h (2, |v). The marginals 
= 1 1 v) of this approximate posterior, together with the data, are used 
to create an augmented input for this deep multilayer neural network, as 
shown in Figure 6. Standard backpropagation of error derivatives can then 
be used to discriminatively fine-tune the model. 9 


9 Note that one can also backpropagate through the "unfolded" mean field, as shown 
in Figure 6, middle panel. In our experiments, however, this did not improve model 
performance. 
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Mean-Field Updates 



Figure 6: (Left) A two-hidden-layer Boltzmann machine. (Right) After learning, 
the DBM model is used to initialize a multilayer neural network. The marginals 
of approximate posterior qihf' = l|v) are used as additional inputs. The net- 
work is fine-tuned by backpropagation. 



The unusual representation of the input is a by-product of converting 
a DBM into a deterministic neural network. In general, the gradient-based 
fine-tuning may choose to ignore Q(h < 2 ) |v), that is, drive the first-layer 
connections W (2> to zero, which will result in a standard neural network. 
Conversely, the network may choose to ignore the input by driving the 
first-layer weights W (i:> to zero and make its predictions based on only the 
approximate posterior. However, the network typically makes use of the 
entire augmented input for making predictions. 

5 Experimental Results 

In our experiments, we used the MNIST and NORB data sets. To speed 
up learning, we subdivided data sets into mini-batches, each containing 
100 cases, and updated the weights after each mini-batch. The number 
of sample particles, used for approximating the model's expected suf- 
ficient statistics, was also set to 100. For the stochastic approximation 
algorithm, we always used five full Gibbs updates of the particles. Each 
model was trained using 300,000 weight updates. The initial learning rate 
was set at 0.005 and was decreased as 10/ (2000 + t), where t is the num- 
ber of updates so far. For discriminative fine-tuning of DBMs, we used 
the method of conjugate gradients on larger mini-batches of 5000 with 
three line searches performed for each mini-batch. (Details of pretrain- 
ing and fine-tuning, along with details of Matlab code that we used for 
learning and fine-tuning deep Boltzmann machines, can be found online at 
http:/ / www.utstat.toronto.edu/~rsalakhu/DBM.html.) 

5.1 MNIST. The MNIST digit data set contains 60,000 training and 
10,000 test images of 10 handwritten digits (0 to 9), with 28x28 pixels. 
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Table 1: Results of Estimating Partition Functions of an RBM, One Fully Con- 
nected BM, and Two DBM Models, Along with the Estimates of the Lower 
Bound on the Average Training and Test Log Probabilities. 



Estimates 

Average Log Probability 

logi 

log ( Z ±a) 

Test 

Train 

RBM 

390.76 

390.56, 390.92 

-86.90 

-84.67 

Flat BM 

198.29 

198.17, 198.40 

-84.67 

-84.35 

2-layer BM 

356.18 

356.06 , 356.29 

-84.62 

-83.61 

3-layer BM 

456.57 

456.34 , 456.75 

-85.10 

-84.49 


Notes: For all BMs we used 20,000 intermediate distributions. Results were averaged over 
100 AIS runs. 


Intermediate intensities between 0 and 255 were treated as probabilities, 
and each time an image was used, we sampled binary values from these 
probabilities independently for each pixel. 

In our first experiment, we trained a fully connected flat BM on the 
MNIST data set. The model had 500 hidden units and 784 visible units. To 
estimate the model's partition function, we used 20,000 /3 k spaced uniformly 
from 0 to 1. Results are shown in Table 1, where for all models, we used 
M = 100 AIS runs (see equation 4.9). The estimate of the lower bound on 
the average test log probability was —84.67 per test case, which is slightly 
better compared to the lower bound of —85.97, achieved by a carefully 
trained two-hidden-layer deep belief network (Salakhutdinov & Murray, 
2008). 

In our second experiment, we trained two deep Boltzmann machines: 
one with two hidden layers (500 and 1000 hidden units) and the other 
with three hidden layers (500, 500, and 1000 hidden units), as shown in 
Figure 8. To estimate the model's partition function, we also used 20,000 
intermediate distributions spaced uniformly from 0 to 1. Table 1 shows 
that the estimates of the lower bound on the average test log probability 
were —84.62 and —85.10 for the two- and three-layer Boltzmann machines, 
respectively. Observe that although the two DBMs contain about 0.9 and 
1.15 million parameters, they do not appear to suffer much from over- 
fitting. The difference between the estimates of the training and test log 
probabilities was about 1 nat. Figure 7 further shows samples generated 
from all three models by randomly initializing all binary states and run- 
ning the Gibbs sampler for 100,000 steps. Certainly all samples look like 
the real handwritten digits. We also emphasize that without greedy pre- 
training, we could not successfully learn good DBM models of MNIST 
digits. 

To estimate how loose the variational bound is, we randomly sampled 
100 test cases, 10 of each class, and ran AIS to estimate the true test log 
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Figure 7: Random samples from the training set, and samples generated from 
three Boltzmann machines by running the Gibbs sampler for 100,000 steps. The 
images shown are the probabilities of the binary visible units given the binary 
states of the hidden units. 
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Figure 8: (Left) The architectures of two deep Boltzmann machines used in 
MNIST experiments. (Right) The architecture of deep Boltzmann machine used 
in NORB experiments. 


probability for the two-layer Boltzmann machine. 10 The estimate of the 
variational bound was —83.35 per test case with an error estimate 
(—83.21. —83.51). The estimate of the true test log probability was —82.86 
with an error estimate (—82.68, 83.12). The difference of about 0.5 nats 
shows that the bound is rather tight. 

For a simple comparison, we also trained several mixture of Bernoullis 
models with 10, 100, 500, 1000, and 2000 components. The corresponding 
average test log probabilities were —168.95, —142.63, —137.64, —133.21, and 


10 Note that computationally, this is equivalent to estimating 100 partition functions, 
as discussed at the end of section 4.1. 
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Table 2: Classification Error Rates on MNIST Test Set When Only a Small Frac- 
tion of Labeled Data Is Available But the Remainder of the 60,000 Training Cases 
Is Available Unlabeled. 



Two-Layer 

DBM 

Nonlinear 

NCA 

Linear 

NCA 

Stack of 
Autoencoders 

KNN 

1% (600) 

4.82% 

8.81% 

19.37% 

8.13% 

13.74% 

5% (3000) 

2.72% 

3.24% 

7.23% 

3.54% 

7.19% 

10% (6000) 

2.46% 

2.58% 

4.89% 

2.98% 

5.87% 

100% (60000) 

0.95% 

1.00% 

2.45% 

1.40% 

3.09% 


—135.78. Compared to DBMs, a mixture of Bernoullis performs very badly. 
The difference of about 50 nats per test case is striking. 

Finally, after discriminative fine-tuning, the two-hidden-layer Boltz- 
mann machine achieves an error rate of 0.95% on the full MNIST test set. 
This is, to our knowledge, the best published result on the permutation- 
invariant version of the MNIST task. 11 The three-layer BM gives a slightly 
worse error rate of 1.01%. The flat BM, on the other hand, gives consid- 
erably worse error rate of 1.27%. This is compared to 1.4% achieved by 
SVMs (Decoste & Scholkopf, 2002), 1.6% achieved by randomly initialized 
backpropagation, 1.2% achieved by the deep belief network, described in 
Hinton et al. (2006) and Hinton and Salakhutdinov (2006), and 0.97% ob- 
tained by using a combination of discriminative and generative fine-tuning 
on the same DBN (Hinton, 2007). 

To test discriminative performance of DBMs when the number of labeled 
examples is small, we randomly sampled 1%, 5%, and 10% of the handwrit- 
ten digits in each class and treated them as labeled data. Table 2 shows that 
after discriminative fine-tuning, a two-hidden-layer BM achieves error rates 
of 4.82%, 2.72%, and 2.46%. Deep Boltzmann machines clearly outperform 
regularized nonlinear NCA (Salakhutdinov & Hinton, 2007), linear NCA 
(Goldberger, Roweis, Hinton, & Salakhutdinov, 2004), a stack of greedily 
pretrained autoencoders (Bengio, Lamblin, Popovici, & Larochelle, 2007), 
and K-nearest neighbors, particularly when the number of labeled exam- 
ples is only 600. We note that similar to DBMs, both regularized nonlinear 
NCA and a stack of autoencoders use an unsupervised pretraining stage 
(trained on all 60,000 unlabeled MNIST images), followed by supervised 
fine-tuning. 

5.2 NORB. Results on MNIST show that deep Boltzmann machines 
can significantly outperform many other models on the well-studied but 


n In the permutation-invariant version, the pixels of every image are subjected to the 
same random permutation, which makes it hard to use prior knowledge about images. 
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Figure 9: Random samples from the training set, and samples generated from 
a three-hidden-layer deep Boltzmann machine by running the Gibbs sampler 
for 10,000 steps. 


relatively simple task of handwritten digit recognition. In this section, we 
present results on NORB, which is a considerably more difficult data set 
than MNIST. NORB (LeCun, Huang, & Bottou, 2004) contains images of 50 
different 3D toy objects with 10 objects in each of five generic classes: cars, 
trucks, planes, animals, and humans. Each object is photographed from 
different viewpoints and under various lighting conditions. The training 
set contains 24,300 stereo image pairs of 25 objects, 5 per class, while the 
test set contains 24,300 stereo pairs of the remaining, different 25 objects. 
The goal is to classify each previously unseen object into its generic class. 
From the training data, 4,300 were set aside for validation. 

Each image has 96x96 pixels with integer grayscale values in the range 
[0,255]. To speed up experiments, we reduced the dimensionality by using 
a foveal representation of each image in a stereo pair. The central 64x64 
portion of an image is kept at its original resolution. The remaining 16 pixel- 
wide-ring around it is compressed by replacing nonoverlapping square 
blocks of pixels in the ring with a single scalar given by the average pixel 
value of a block. We split the ring into four smaller ones: the outermost 
ring consists of 8x8 blocks, followed by a ring of 4x4 blocks, and finally 
two innermost rings of 2x2 blocks. The resulting dimensionality of each 
training vector, representing a stereo pair, was 2 x 4488 = 8976. A random 
sample from the training data used in our experiments is shown in Figure 9 
(left). 

To model raw pixel data, we use an RBM with gaussian visible and 
binary hidden units. Gaussian RBMs have been previously successfully 
applied for modeling greyscale images, such as images of faces (Hinton & 
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Salakhutdinov, 2006). However, learning an RBM with gaussian units can 
be slow, particularly when the input dimensionality is quite large. Here we 
follow the approach of Nair and Hinton (2009) by first learning a gaussian 
RBM and then treating the activities of its hidden layer as "preprocessed" 
data. Effectively, the learned low-level RBM acts as a preprocessor that 
converts grayscale pixels into a binary representation, which we then use 
for learning a deep Boltzmann machine. 

The number of hidden units for the preprocessing RBM was set to 4000, 
and the model was trained using contrastive divergence learning for 500 
epochs. We then trained a two-hidden-layer DBM with each layer contain- 
ing 4000 hidden units, as shown in Figure 8 (right). Note that the entire 
model was trained in a completely unsupervised way. After the subsequent 
discriminative fine-tuning, the "unrolled" DBM achieves a misclassification 
error rate of 10.8% on the full test set. This is compared to 11.6% achieved 
by SVMs (Bengio & LeCun, 2007), 22.5% achieved by logistic regression, 
and 18.4% achieved by the K-nearest neighbors (LeCun et al., 2004). To 
show that DBMs can benefit from additional unlabeled training data, we 
augmented the training data with additional unlabeled data by applying 
simple pixel translations, creating 1,166,400 training instances. 12 After learn- 
ing a good generative model, the discriminative fine-tuning (using only 
the 24,300 labeled training examples without any translation) reduces the 
misclassification error to 7.2%. Figure 9 shows samples generated from the 
model by running prolonged Gibbs sampling. Note that the model was able 
to capture a lot of regularities in this high-dimensional, richly structured 
data, including different object classes, various viewpoints, and lighting 
conditions. 

Finally, we tested the ability of the DBM to perform an image in-painting 
task. To this end, we randomly selected 10 objects from the test set and sim- 
ulated the occlusion by zeroing out the left half of the image (see Figure 10). 
We emphasize that the test objects are different from the training objects 
(i.e., the model never sees images of "cowboy," but it sees other images 
belonging to the "person" category). We next sampled the "missing" pixels 
conditioned on the nonoccluded pixels of the image using 1000 Gibbs up- 
dates. Figure 10 (bottom) shows that the model was able to coherently infer 
occluded parts of the test images. In particular, observe that even though the 
model never sees an image of the cowboy, it correctly infers that it should 
have two legs and two arms. 

Surprisingly, even though the deep Boltzmann machine contains about 
68 million parameters, it significantly outperforms many of the competing 
models. Clearly unsupervised learning helps generalization because it en- 
sures that most of the information in the model parameters comes from 
modeling the input data. The very limited information in the labels is used 


12 We thank Vinod Nair for sharing his code for blurring and translating NORB images. 


2000 


R. Salakhutdinov and G. Hinton 


* 


\ 

L 

£ 


% 

> 


* 


1 

■ 

w 

0 


** 


M' 

m 

L 

■ 

1 

u 

“fca 

> 


\ 




V 

y 

— 

> 


Figure 10: Performance of the three-hidden-layer DBM on the image in- 
painting task. (Top) Ten objects randomly sampled from the test set. 
(Middle) Partially occluded input images. (Bottom) Inferred images were gene- 
rated by running a Gibbs sampler for 1000 steps. 


only to slightly adjust the layers of features already discovered by the deep 
Boltzmann machine. 

6 Discussion 

A major difference between DBNs and DBMs is that the procedure for 
adding an extra layer to a DBN replaces the whole prior over the previous 
top layer, whereas the procedure for adding an extra layer to a DBM replaces 
only half of the prior. So in a DBM, the weights of the bottom-level RBM 
end up doing much more of the work than in a DBN where the weights are 
used only to define p(v|h (1) ; W 1 ) (in the composite generative model). This 
suggests that adding layers to a DBM will give diminishing improvements 
in the variational bound much more quickly than adding layers to a DBN. 
There is, however, a simple way to pretrain a DBM so that more of the 
modeling work is left to the higher layers. 

Suppose we train an RBM with one set of hidden units, and four sets of 
visible units, and we constrain the four weight matrices (and visible biases) 
to be identical. 13 Then we use the hidden activities as data to train an RBM 
with one set of visible units and four sets of hidden units, again constrained 
to have identical weight matrices. Now we can combine the two RBMs into 
a DBM with two hidden layers by using one copy of the weight matrix from 
the first RBM and three times one of the copies of the weight matrix from 
the second RBM. In this DBM, three-fourths of the first RBM's prior over 
the first hidden layer has been replaced by the prior defined by the second 
RBM. It remains to be seen whether this makes DBMs work better. It is also 
not obvious how this idea can be applied to intermediate hidden layers. 


13 As before, we use mean-field reconstructions of the four sets of visible units to avoid 
modeling the fact that all four sets have the same states in the data. 
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In this article, we have focused on Boltzmann machines with binary 
units. The learning methods we have described can be extended to learn 
deep Boltzmann machines built with RBM modules that contain real-valued 
(Marks & Movellan, 2001), count (Salakhutdinov & Hinton, 2009b), or tab- 
ular data provided the distributions are in the exponential family (Welling, 
Rosen-Zvi, & Hinton, 2005). However, it often requires additional insights 
to get the basic RBM learning module to work well with nonbinary units. 
For example, it ought to be possible to learn the variance of the noise model 
of the visible units in a Gaussian-Bernoulli RBM, but this is typically very 
difficult for reasons explained in Hinton (2010). For modeling the NORB 
data, we used fixed variances of 1, which is clearly much too big for data 
that have been normalized so that the pixels have a variance of 1. Recent 
work shows that gaussian visible units work much better with rectified lin- 
ear hidden units (Nair & Hinton, 2010) and using this type of hidden unit 
it is straightforward to learn the variance of the noise model of each visible 
unit. 

7 Summary 

We presented a novel combination of variational and Markov chain Monte 
Carlo algorithms for training Boltzmann machines. When applied to pre- 
trained deep Boltzmann machines with several hidden layers and millions 
of weights, this combination is a very effective way to learn good genera- 
tive models. We demonstrated the performance of the algorithm using the 
MNIST handwritten digits and the NORB stereo images of 3D objects with 
highly variable viewpoint and lighting. 

A simple variational approximation works well for estimating the data- 
dependent statistics because learning based on these estimates encourages 
the true posterior distributions over the hidden variables to be close to 
their variational approximations. Persistent Markov chains work well for 
estimating the data-independent statistics because learning based on these 
estimates encourages the persistent chains to explore the state space much 
more rapidly than would be predicted by their mixing rates. 

Pretraining a stack of RBMs using contrastive divergence can be used 
to initialize the weights of a deep Boltzmann machine to sensible values. 
The RBMs can then be composed to form a deep Boltzmann machine. The 
pretraining ensures that the variational inference can be initialized sensibly 
by a single bottom-up pass from the data vector using twice the bottom-up 
weights to compensate for the lack of top-down input on the initial pass. 

We further showed how annealed importance sampling, along with vari- 
ational inference, can be used to estimate a variational lower bound on 
the log probability that a deep Boltzmann machine assigns to test data. 
This allowed us to directly assess the performance of deep Boltzmann ma- 
chines as generative models of data. Finally, we showed how to use a deep 
Boltzmann machine to initialize the weights of a feedforward neural 
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network that can then be discriminatively fine-tuned. These networks give 
excellent discriminative performance, especially when there are very few 
labeled training data but a large supply of unlabeled data. 

Appendix 

For clarity of presentation, let p be a shorthand notation for the distribu- 
tion P(h (1) ; W (1) ), defined in equation 3.12, and let p h denote an individual 
term from this distribution. Similarly, we will let r h denote an individual 
term from P(h (1) ; W (2) ), q h denote a term from Q(h ( 1 ) |v n ; W (1) ), and m h de- 
note a term from the geometric mean of the two probability distributions 
P(h (1) ; W fl) ), P(h (1) ; W (2) ), 


which is the renormalized pairwise product of the square roots of the two 
probabilities for each event (see equation 3.14). We note that the normalizing 
constant Z is equal to 1 if the two distributions p and r are identical, and it 
is less than 1 otherwise. 

For each case n, given 



(A.l) 


h 


KLOfU r) < KL(q\\p), 


(A. 2) 


we want to show that KL(q\\m ) < KL(q\\p): 


KL(q\\m) = J^q h log^ 


h 




(follows from equation A.2) 


< J2 h lo § % - J2 1h lo § Vh + lo g z 


h 


h 


J2 h lo § % - J2 1h lo § Vh (since Z < 1) 


h 


h 


<KL(q\\p). 
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The above derivation is a special case of a more general result, provided 
in Hinton (2002). It states that the Kullback-Leibler divergence between the 
geometric mean of a set of probability distributions P : and the distribution 
Q is smaller than the average of the Kullback-Leibler divergences of the 
individual distributions: 


where w i are nonnegative and sum to 1. Note that the normalizing constant 

Z is 1 if all of the individual distributions are identical. Otherwise, Z < 1, 

and the difference between the two sides of the above equation is log(l/.Z). 
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